Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor the Tian Approximation in the Reactive Fluid Transport Model #6161

Open
wants to merge 1 commit into
base: main
Choose a base branch
from

Conversation

danieldouglas92
Copy link
Contributor

I was looking through the Reactive fluid transport model code and saw how the katz 2003 model is in it's own reaction model. I think I recall at the hackathon that the ultimate goal for this material model was to expand it by creating new reaction models, so I went and refactored the tian approximation into it's own reaction.

Copy link
Member

@gassmoeller gassmoeller left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for the delayed review. This looks generally good, I just have a few comments. Please add a changelog entry that mentions this plugin is now easier to include in material models.

Comment on lines +21 to +22
#ifndef _aspect_material_reaction_melt_tian2019_solubility_h
#define _aspect_material_reaction_melt_tian2019_solubility_h
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know you took this from the other reaction model, but that one sets a bad example here. These lines should correspond as closely as possible to the folder structure and file name. This is not technically necessary, but a convention we follow to make them easy to remember:

Suggested change
#ifndef _aspect_material_reaction_melt_tian2019_solubility_h
#define _aspect_material_reaction_melt_tian2019_solubility_h
#ifndef _aspect_material_model_reaction_model_tian2019_solubility_h
#define _aspect_material_model_reaction_model_tian2019_solubility_h


template <int dim>
Tian2019Solubility<dim>::Tian2019Solubility()
= default;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you make the constructor the default constructor then you usually do not need to implement it at all. Just remove it from here and from the header file.

Comment on lines +57 to +68
/**
* Declare the parameters this function takes through input files.
*/
static
void
declare_parameters (ParameterHandler &prm);

/**
* Read the parameters from the parameter file.
*/
void
parse_parameters (ParameterHandler &prm);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please put these functions at the end of the list of public functions here and in the source file. This is just a convention we (try to) follow, but it makes navigating the source files easier.
A typical class structure in ASPECT looks like this (with many functions being optional):


class ClassName
{
  constructor();

  initialize();

  all_other_minor_functions();

  main_function(in this case melt_fraction);

  declare_parameters();

  parse_parameters();

  private:

  ... all member variables and private functions
}



/**
* Compute the free fluid fraction that can be present in the material based on the
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can be present or is present? I think you mean is present because can be present would be 100%?

/**
* Compute the free fluid fraction that can be present in the material based on the
* fluid content of the material and the fluid solubility for the given input conditions.
* @p in and @p melt_fractions need to have the same size.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

there is no melt_fractions variable

Comment on lines +90 to +91
* @param q unsigned int from 0-3 indexing which rock phase the equilbrium
* bound water content is being calculated for
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, is that really the meaning of q? We usually use q as an index for the quadrature point to index into the vectors that are contained inside in. Please clarify and if you use this variable for something else, maybe find a different name unless it has to be q for other reasons.

Comment on lines +97 to +99
/**
* Parameters for the solubility model of Tian et al., 2019.
*/
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this comment here is redundant now. All of these parameters are from Tian et al since this is the Tian et al plugin.

Comment on lines +119 to +133
std::vector<double> LR_peridotite_poly_coeffs {-19.0609, 168.983, -630.032, 1281.84, -1543.14, 1111.88, -459.142, 95.4143, 1.97246};
std::vector<double> csat_peridotite_poly_coeffs {0.00115628, 2.42179};
std::vector<double> Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603};

std::vector<double> LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519};
std::vector<double> csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732};
std::vector<double> Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517};

std::vector<double> LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365};
std::vector<double> csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588};
std::vector<double> Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049};

std::vector<double> LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459};
std::vector<double> csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867};
std::vector<double> Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All of these numbers are fixed and will never change during the program runtime. You can tell this to the compiler, which will then optimize the code further and also prevent you from accidentally changing the values by changing their type from std::vector<double> to constexpr std::array<double,N> where N is the number of values in this array. Initialization and access to these arrays can remain unchanged.

* and 50 GPa for sediment. These cutoff pressures were determined by extending the pressure range in Tian et al. (2019)
* and observing where the maximum allowed water contents jump towards infinite values.
*/
const std::vector<double> pressure_cutoffs {10, 26, 16, 50};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here, make this a constexpr std::array

Comment on lines +367 to +372
prm.enter_subsection("Tian 2019 model");
{
tian2019_model.initialize_simulator (this->get_simulator());
tian2019_model.parse_parameters(prm);
}
prm.leave_subsection();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So you put the parameters of the Tian model in their own subsection now. However, this means the parameter files of everyone using the reactive fluid flow model with the 3.0 version will need to be updated. Please include a changelog entry that mentions this. Additionally, you could extend our update scripts in contrib/utilities/update_scripts/prm_files/ to automatically move this parameter. The closest operation that does something like this is probably this function. However, the update process looks a bit cryptic, I leave this up to you to decide if this is important enough.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants